"""McCrackin's method for solving psi-delta pairs to give refractive index and film thickness"""
"""This file is part of PyEllips.

    PyEllips is free software: you can redistribute it and/or modify
    it under the terms of the GNU General Public License as published by
    the Free Software Foundation, either version 3 of the License, or
    (at your option) any later version.

    PyEllips is distributed in the hope that it will be useful,
    but WITHOUT ANY WARRANTY; without even the implied warranty of
    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
    GNU General Public License for more details.

    You should have received a copy of the GNU General Public License
    along with PyEllips.  If not, see <http://www.gnu.org/licenses/>."""
#import fresnel as fr
import cmath
import numpy as np
import scipy as sp

from Tkinter import *
import tkMessageBox
class Ellipsometer:
	"Class to hold ellipsometry parameters-angle, wavelength, substrate RI etc"
	def __init__(self,angle,wavel,ns,ks):
		self.angle=angle
		self.wavel=wavel
		self.subsri=complex(ns,ks)
		self.angler=np.pi*angle/180.0
		#tkMessageBox.showinfo("Text", "angle =%f  angler=%f" % (angle,self.angler)) 
		
		
	def find_nd(self,nstart,nend,psi,delta):
		#find best combination of thickness and refractive index in range nstart to nend
		bestd=self.calcnd(psi,delta,nstart)
		#print bestd.imag
		bestn=nstart
		for i in np.arange(nstart+0.01,nend,0.005):
			#print i
			test=self.calcnd(psi,delta,i)
			if np.fabs(bestd.imag)>np.fabs(test.imag):
				#print np.fabs(bestd.imag),np.fabs(test.imag)
				bestd=test
				bestn=i
		# have been through the values and have optimum value in best
		return (bestn,bestd)
		
	def calcnd(self,psi,delta,testn):
	#returns the thickness value that the McCrackin method 
	#gives for the input psi delta and testn."""
		#setup refractive indices of layers
		n1=complex(1.0,0.0) #air
		n2=complex(testn,0.) # film value being tried
		#need to make these user definable
		#n3=self.subsri.get()
		n3=self.subsri		
		#convert psi,delta to radians
		psir=np.pi*psi/180.
		deltar=np.pi*delta/180.
		tanp=cmath.tan(complex(psir,0.))
		#calculate the angles in the layers
		cos1=complex(np.cos(self.angler),0.)
		cos2=self.cos_snell(n1,n2,self.angler)
		cos3=self.cos_snell(n1,n3,self.angler)
		
# now calculate the fresnel terms
	# rp
		rp12=self.rp(n1,n2,cos1,cos2)
		rp23=self.rp(n2,n3,cos2,cos3);
	
	#rs
		rs12=self.rs(n1,n2,cos1,cos2);
		rs23=self.rs(n2,n3,cos2,cos3);
		
	# quadratic terms
	
		A=tanp*cmath.exp(deltar*1j)*rp12*rp23*rs23-rp23*rs12*rs23
	

		B=tanp*cmath.exp(deltar*1j)*(rs23+rp12*rp23*rs12)-(rp12*rs12*rs23+rp23)
		C=tanp*cmath.exp(deltar*1j)*rs12 - rp12
		#print A,B,C
	# solutions to quadratic equation
		X1=(-B+cmath.sqrt((B*B)-(4.*A*C)))/(2.*A)
		   
		X2=(-B-cmath.sqrt((B*B)-(4.*A*C)))/(2.*A)
 		# convert these values into thicknesses 
		d2a=((self.wavel*1j)/(4.*np.pi*n2*cos2))*cmath.log(X1);

		d2b=((self.wavel*1j)/(4.*np.pi*n2*cos2))*cmath.log(X2);
	
	# return the value with the smallest imaginary part
		if np.fabs(d2a.imag)<np.fabs(d2b.imag):
			return d2a
		else:
			return d2b
			
	def calcpsiDelta(self,nfilm,dfilm):
		#Returns psi and delta values for given n,d
		#setup refractive indices of layers
		n1=complex(1.0,0.0) #air
		n2=complex(nfilm,0.0) # film value being tried
		#need to make these user definable
		#n3=self.subsri.get()
		n3=self.subsri		
		
		#calculate the angles in the layers
		cos1=complex(np.cos(self.angler),0.)
		cos2=self.cos_snell(n1,n2,self.angler)
		cos3=self.cos_snell(n1,n3,self.angler)
		#print cos1,cos2,cos3
# now calculate the fresnel terms
	# rp
		rp12=self.rp(n1,n2,cos1,cos2)
		rp23=self.rp(n2,n3,cos2,cos3);
		#print
		#print rp12,rp23
	
	#rs
		rs12=self.rs(n1,n2,cos1,cos2);
		rs23=self.rs(n2,n3,cos2,cos3);
		#print rs12,rs23
		#expD=cmath.exp(-4.j*np.pi*nfilm*cmath.cos(self.angler)*dfilm/self.wavel)
		expD=cmath.exp(-4.j*np.pi*nfilm*cos2*dfilm/self.wavel)
		
		Rp=(rp12+rp23*expD)/(1+rp12*rp23*expD)
		Rs=(rs12+rs23*expD)/(1+rs12*rs23*expD)
		rho=Rp/Rs
		#print Rp,Rs
		#print rho
		#print np.abs(rho)
		#print abs(rho)
		psi=np.arctan(abs(rho))
		delta=cmath.phase(rho)
		#delta=np.atan2(rho.imag,rho.real)
		return (self.radianstodegrees(psi),self.radianstodegrees(delta))

	def degreestoradians(self,degrees):
		return (np.pi*degrees/180.0)
	
	def radianstodegrees(self,radians):
		return (180.0*radians/np.pi)
	
	def cos_snell(self,n1,n2,angle):
#uses a combination of Snell's law and 1=cos^2+sin^2 to calculate the cosine of the incident angle in the given media
		return(cmath.sqrt(1.0-(cmath.sin(angle)*(n1/n2))**2.));


	def rp(self, n1, n2,  cangle1,  cangle2):

	# calculates the fresnel term rp
	# the cangle variables are the cosines of the angles.
		return ((n2*cangle1-n1*cangle2)/(n2*cangle1+n1*cangle2))
	
	def rs(self,n1,n2,cangle1,cangle2):

	# calculates the fresnel term rs
	# the cangle variables are the cosines of the angles.
		return ((n1*cangle1-n2*cangle2)/(n1*cangle1+n2*cangle2))
		
 
		
		

